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1.  SUMMARY 


The  ability  to  pre-1  >t  the  coupled  effects  of  complex  transport  phenomena  with  de¬ 
tailed  chemical  kinetic;,  is  critic?  in  an  understanding  of  turbulent  reacting  flows,  in  im¬ 
proving  engine  efficiency,  and  in  the  study  of  pollutant  formation.  Since  three-dimensional 
models  combining  both  fluid  dynamical  effects  with  finite  rate  chemistry  are  as  yet  com¬ 
putationally  infeasible,  the  modeling  of  chemically  reacting  flows  generally  proceeds  along 
two  independent  paths.  In  one  case  chemistry  is  given  priority  over  fluid  mechanical  effects 
and  these  models  are  used  to  assess  the  important  elementary  reaction  paths,  for  exam¬ 
ple,  in  hydrocarbon  fuels.  In  the  other  case,  multidimensional  fluid  dynamical  effects  are 
emphasized  with  chemistry  receiving  little  or  no  priority. 

The  goal  in  reacting  flow  computations  is  to  be  able  to  combine  the  effects  of  detailed 
chemistry  with  complex  fluid  mechanics.  During  the  past  six  months,  our  Phase  I  SBIR 
concentrated  in  this  direction.  We  investigated  several  areas  in  which  the  modeling  of 
chemically  reacting  flows  could  be  improved: 

1.  We  incorporated  detailed  kinetics  with  complex  transport  in  the  solution  of  one¬ 
dimensional  diffusion  flames. 

2.  We  investigated  flame  sheet  initialization  procedures. 

3.  We  considered  more  implicit  solution  algorithms  for  the  highly  nonlinear  discrete  equa¬ 
tions. 


As  a  result  of  the  limited  computer  resources  available  under  Phase  I  combined  with 
the  computational  complexity  of  calculating  chemically  reacting  flows,  we  devoted  a  large 
part  of  our  effort  to  the  solution  of  counterflow  diffusion  flames.  Unlike  models  in  which 
the  chemistry  and  transport  phenomena  are  approximated  by  simplified  relations,  we  in¬ 
corporated  both  detailed  finite  rate  chemistry  and  complex  kinetic  theory  transport  effects 
in  our  model.  We  then  investigated  the  use  of  a  flame  sheet  model  as  a  means  of  providing 
improved  starting  estimates  for  the  finite  rate  counterflow  problem.  Solution  of  the  flame 
sheet  and  finite  rate  equations  was  accomplished  with  a  combination  of  time  integration 
and  Newton’s  method.  Based  upon  the  success  of  our  study  in  one  dimension,  we  directed 
our  efforts  to  two-dimensional  problems.  In  particular,  we  generalized  our  one-dimensional 
flame  sheet  model  to  two  dimensions  and,  in  parallel  with  this  effort,  we  investigated  ways 
of  improving  the  computational  efficiency  of  the  primitive  variable  implementation  of  the 
TEACH  code.  The  results  of  our  Phase  I  research  clearly  indicate  the  feasibility  of  our 
methodology  in  solving  reacting  systems  with  detailed  transport  and  complex  chemistry. 
This  methodology  has  potential  applications  to  problems  in  turbulent  reacting  flows,  en¬ 
gine  efficiency,  commercial  power  generation  units,  and  pollutant  formation. 


Dist 


2.  COUNTERFLOW  DIFFUSION  FLAMES 


The  flame  type  of  most  practical  combvistion  devices  is  the  diffusion  flame.  These 
flames  are  important  in  the  interaction  of  heat  and  mass  transfer  with  chemical  reactions  in 
ram  jets,  jet  turbines  and  commercial  burners.  Three-dimensional  models  that  couple  the 
effects  of  fluid  flow  with  detailed  chemical  reactions  are  as  yet  computationally  infeasible. 
We  can,  however,  obtain  important  information  in  practical  systems  by  considering  two- 
dimensional  configurations.  The  counterflow  diffusion  flame  is  one  such  configuration. 

Experimentally  these  flames  can  be  produced  when  a  reaction  zone  is  stabilized  near 
the  stagnation  point  of  two  infinitely  wide  coaxial  concentric  jets  [1]  (see  Figure  I).  Fuel 
is  emitted  from  one  jet  and  oxidizer  (air)  from  the  other.  Combustion  occurs  within  a 
thin  flame  zone  near  the  stagnation  point  where  the  fuel  and  the  oxidizer  are  in  stoi¬ 
chiometric  proportion.  Although  the  flow  in  the  double  jet  experimental  configuration  is 
two-dimensional,  the  mathematical  model  can  be  reduced  to  the  solution  of  a  system  of 
coupled  nonlinear  two-point  boundary  value  problems  along  the  stagnation  point  stream¬ 
line.  In  this  way  we  can  investigate  the  effect  of  detailed  chemical  kinetics  with  complex 
transport  while  still  having  a  computationally  feasible  problem. 


2.1  Problem  Formulation 

Our  model  for  counterflow  diffusion  flames  assumes  the  flow  to  be  laminar,  stagnation 
point  flow.  Hence,  the  governing  boundary  layer  equations  for  mass,  momentum,  chemical 
species  and  energy  can  be  written  in  the  form 

d(puxa)  d(pvxa)  n  . 

~i^+Ai>r=0-  <2-‘> 

du  du  dp  d  (  du\  .  . 

pUTx+pVd-y  +  Tx  =  d-y\!ld-y)'  (22) 

dYk  dYk  d  , 

Pu~gJ  +  pv~Q^  +  fa}  ( pYkVky )  -  Wkwk  =  0,  k  =  1,2,...,  K<  (2.3) 

dT  dT  d  (d  T\  A  dT  A 

^ a7  +  ~  Tv  +  L  +  E  'h'w“h'  =  °-  W 


where  a  represents  a  geometric  factor  (a  =  0  for  cartesian  coordinates  and  o  =  I  for 
cylindrical  coordinates).  The  system  is  closed  with  the  ideal  gas  law, 


In  these  equations  x  and  //  denote  independent  spatial  coordinates  in  the  tangential 
and  transverse  directions,  respectively;  T.  the  temperature;  Yk ,  the  mass  fraction  of  the 


k1^  species;  p,  the  pressure;  u  and  v  the  tangential  and  the  transverse  components  of  tin* 
velocity,  respectively;  p,  the  mass  density;  Wk,  the  molecular  weight  of  the  k^1  species;  TT. 
the  mean  molecular  weight  of  the  mixture;  R ,  the  universal  gas  constant;  A,  the  thermal 
conductivity  of  the  mixture;  c,,,  the  constant  pressure  heat  capacity  of  the  mixture;  cPk. 
the  constant  pressure  heat  capacity  of  the  k^  species;  wk,  the  molar  rate  of  production  of 
the  A:1*1  species  per  unit  volume;  hk,  the  specific  enthalpy  of  the  k ^  species;  p  the  viscosity- 
of  the  mixture  and  Vky  is  the  diffusion  velocity  of  the  k ^  species  in  the  y  direction.  In 
both  configurations  the  free  stream  (tangential)  velocity  at  the  edge  of  the  boundary  layer 
is  given  by  =  ax  where  a  is  the  strain  rate. 


Upon  introducing  the  notation 


/'=  — , 
uoo 

V  =  pv. 


the  boundary  layer  equations  can  be  transformed  into  a  system  of  ordinary  differential 
equations  valid  along  the  stagnation-point  streamline  x  =  0.  For  a  system  in  rectangular 
or  cylindrical  coordinates,  we  have 


dV 

—  +  a(l  +  a)pf  =  0. 

dy 

^  ("^)  -v'^r-t-“(eoo-^(/')2)  =  o. 

~~{pYkYk)  -Y^  +  wkWk  =  0,  k  —  1,2 . A , 

dy  dy 


d  (  dT\  ..dT  A  ....  dT  o 

Ty  \~dy  )  "  CpV  d^  ~  L.M*' kTy  "  £  WkW'hk  =  °‘ 


(2.10) 


(2.11) 


The  boundary  conditions  for  the  double-jet  configuration  at  y  =  —  oo  are  given  by 


V  —  V 

v  —  V  —  OO* 


(2.12) 


and  at  y  =  x>  by 


n  =  n_,  k  — 1,2 . a 

T  =  H*, 


f  =  1, 

Yk  =  Yu,  k  —  1.2 . A. 
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(2.18) 


1 


T  =  Too. 

The  mass  flux,  temperature,  and  species  mass  fractions  (l'_  T_oo,  V*_ x, )  at  the  fuel  jet 
are  specified  quantities  as  are  the  temperature  and  species  mass  fractions  {Too  and  Yk,^ ) 
at  the  oxidizer  jet. 

2.2  Flame  Sheet  Model 

Despite  the  outwardly  simple  form  of  the  counterflow  equations,  the  determination  of 
a  "good”  initial  solution  estimate  can  be  difficult.  The  difficulty  is  due  to  the  exponential 
dependence  of  the  chemistry  terms  on  the  temperature  and  to  the  nonlinear  coupling  of 
the  fluid  and  the  thermochemistry  solution  fields.  In  previous  wrork,  we  focused  our  efforts 
on  the  solution  of  adiabatic  and  nonadiabatic  premixed  laminar  flames  by  adaptive  finite 
difference  methods  [2-3) .  In  these  problems  cubic  polynomials  and  Gaussian  shaped  profiles 
were  used  as  starting  estimates  for  the  major  and  minor  species  on  an  initial  coarse  grid. 
These  approximations  were  often  sufficient  to  bring  the  starting  estimates  into  the  domain 
of  convergence  of  Newton’s  method. 

In  adiabatic  and  nonadiabatic  premixed  laminar  flame  problems  the  conservation  of 
mass  and  momentum  reduces  to  the  specification  of  a  constant  mass  flow  rate  and  a 
constant  thermodynamic  pressure  [2-3].  Hence,  thermochemical  considerations  play  a  more 
important  role  in  these  problems  than  do  fluid  dynamical  aspects.  This  is  not  the  case 
in  counterflow  diffusion  flames.  In  particular,  there  is  a  strong  coupling  between  the 
fluid  dynamic  and  the  thermochemistry  solution  fields  in  these  flames.  Wo  have  found 
that,  although  the  solution  procedure  used  in  premixed  laminar  flame  problems  can  work 
in  selected  counterflow  cases,  it  does  not  provide  a  sufficiently  robust  or  efficient  starting 
estimate  from  which  Newton’s  method  will  converge.  In  addition,  the  relaxation  to  steady- 
state  (or  at  least  until  the  solution  is  within  the  convergence  domain  of  Newton's  method)  is 
very  slow.  The  importance  of  these  flames  in  turbulence  modeling  and  in  the  determination 
of  chemically  controlled  extinction  limits,  however,  necessitates  the  development  of  an 
efficient  starting  procedure.  We  couple  the  appropriate  equations  of  mass  and  momentum 
with  a  Shvab-Zeldovieh  equation  to  provide  flame  sheet  starting  estimates  for  the  mass 
flux  in  the  transverse  direction,  the  similarity  function,  the  temperature,  and  the  stable 
major  species  in  the  flame. 

Our  starting  point  is  the  assumption  that  the  fuel  and  the  oxidizer  obey  a  single  overall 
irreversible  reaction  of  the  type 

Fuel  (/•')  +  Oxidizer  (O)  —  Products  (P),  (2.19) 

in  the  presence  of  an  inert  gas  (N).  We  have 

ufF  +  VuO  —  uPP.  (2.20) 
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where  vp,vo  and  vP  are  the  stoichiometric  coefficients  of  the  fuel,  the  oxidizer  and  the 
product,  respectively  In  addition,  we  neglect  thermal  diffusion  and  assume  that  the  ordi¬ 
nary  mass  diffusion  v<  'ooities  can  be  written  in  terms  of  Fick’s  law,  that  is 


DkdYk 

lit  =  ~— — 5 — ,  k  =  1,2,  ...,A, 
**  ay 


(2  21) 


where  D*  is  the  diffusion  coefficient  of  the  species  with  respect  to  the  mixture.  We 
also  take  cp  —  cPk  to  be  constant.  With  these  approximations  we  can  write 


dV 

dy 


+  a(l  -I-  a)pf'  =  0, 


d_ 

dy 


+  a(?°°  ~  =  °- 

±  (  n  dYF\  iYF 
dy  \  F  dy  )  dy 


Wpi/pw  —  0, 


d  (  dY0\  dY0 

—  { pD0—j—  )  ~  L  — r - W0v0u>  =  0, 


dy 

d_ 

dy 


dy 


dy 


pDp 


dYP\  ,rdYP 

-j—  ~  V  — - 1-  WpUpw  =  0, 

dy  )  dy 


(2.22) 

(2.23) 

(2.24) 

(2.25) 

(2.26) 

(2.27) 

(2.28) 

(2.29) 

Up  Vo  Vp 

is  the  rate  of  progress  of  the  reaction  and  where  we  have  made  use  of  the  fact  that 

£f„nvi  =  o. 

If  we  introduce  the  heat  release  per  unit  mass  of  the  fuel  Q  where 


±(  n  dYN\  dYN 
dy  X  N  dy  )  dy 


d  f  A  dT  \  .  dT  (WpVphp  +  WoVoho  ~  WpVphp) 

— - _  i  _  v  -| - -w  =  0, 

dy  \Cpuy )  dy  cp 


where 


y,  _  WF  —  w°  _  Wp 


ri  —  L  ,  WOuO  u  WpVp 
Q  —  hp  +  777— — ho  -  — - hp. 


W/  pVp 

and  if  we  assume  that  the  Lewis  numbers 

A 


W'  pVp 


Ler  = 


Leo  = 


pDFcp 

A 

pDpcp 


.Le0  = 

,  Lcn  = 


pD0cp  ’ 

A 

pDy  Cp 


(2.30) 


(2.31a) 


(2.315 


(2.32) 


are  all  equal  to  one,  then  each  of  the  Shvab-Zeldovich  variables 

Zf  =  Yf-  5>„  +  |(T  -  r„). 


Zo  =  Yo  T„), 

(2.33) 

(2.34) 

II 

£ 

1 

£ 

8 

(2.35) 

satisfies  the  following  differential  equation 

d(  dZk  \  dZk  .. 

—  IpDk-j—  )  -  V =  0,  k  =  F,0,P,N , 
dy  V  dy  )  dy 

(2.36) 

with 

Zk (  00 )  =  Zk_x, 

(2.37) 

Zk(oo)  =  0, 

(2.38) 

for  the  double-jet  problem  where  Zk_OQ  is  constant.  As  a  result,  all  of  the  Zk  are  propor¬ 
tional  to  each  other  and  to  the  conserved  scalar  S  which  satisfies 


5(— oo)  =  1, 

S(oo)  =  0, 

where  D  is  a  diffusion  coefficient. 

From  (2.36)  and  (2.39)  we  can  write 

Zk  =  Zk_xS(y ),  k  =  F,0,  P,  N .  (2.42) 


(2.40) 

(2.41) 


Equation  (2.39)  can  be  coupled  with  equations  (2.22)  and  (2.23)  to  obtain  profiles  for  V.  f 
and  5.  To  complete  the  specification  of  the  starting  estimate,  we  must  be  able  to  recover 
the  temperature  and  the  major  species  profiles  from  the  conserved  scalar.  If  we  utilize  the 
result  in  (2.42)  along  with  the  fact  that  in  the  flame  sheet  model  fuel  and  oxidizer  cannot 
co-exist,  we  can  derive  relations  for  the  temperature  and  major  species  on  either  side  of 
the  reaction  zone.  On  the  fuel  side,  we  have 

r=r_„5+[r„  +  r0„^2^](i-5),  (2.43) 

CpWoL'O. 


n  =  yfs  +  r0c 


IV  fUjt 
W0v0 


(2.44) 
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(2.46) 


and 


YN  =  Y„ao(l-S)  +  YN_ooS. 


On  the  oxidizer  side,  we  have 

T  =  Too(l  -  5)  + 

=0, 


Q 

-Yr  +  T_c 


LCP 


5, 


and 


»  = 

F/v  =  Fv00(l-5)  +  y)v_oo5. 


(2.47) 

(2.48) 

(2.49) 

(2.50) 

(2.51) 

(2.52) 


The  diffusion  coefficient  and  the  viscosity  can  be  determined  by  specifying  a  reference 
Prandtl  number  and  a  transport  relation  for  the  viscosity.  In  particular,  we  let  (Pr)re/  = 
0.75  (air)  and 


(2.53) 


where  r  =  0.7,  T0  =  298  K  and  /i0  =  1.85  x  10-4  gm/cm-sec  is  again  a  reference  value  for 
air.  The  scaled  heat  release  parameter  Q/cp  =  A T  can  be  determined  from  the  heat  of 
combustion  of  the  system  and  a  representative  heat  capacity. 


2.3  Method  of  Solution 

Solution  of  the  governing  equations  proceeds  with  an  adaptive  nonlinear  boundary 
value  method.  Our  goal  is  to  obtain  a  discrete  solution  of  the  governing  equations  on  the 
mesh  X.  where 

X  =  {-L  =  y0  <  tji  <  ■  ■ .  <  ym  =  L}.  (2.54) 

With  the  continuous  differential  operators  replaced  by  finite  difference  expressions,  we 
convert  the  problem  of  finding  an  analytic  solution  of  the  governing  equations  to  one  of 
finding  an  approximation  to  this  solution  at  each  point  of  the  mesh  X.  We  seek  the  solution 
of  the  nonlinear  system  of  difference  equations 


For  an  initial  solution  estimate  IJ°  that  is  sufficiently  “close”  to  U£,  the  system  of  equations 
in  (2.55)  can  be  solved  by  Newton's  method.  We  write 

J{Uk)  (Uk+l  -  Uk)  =  -\kF(Uk),  k  =  0,1 .  (2.56) 

where  Uk  denotes  the  k^1  solution  iterate.  A*  the  k ^  damping  parameter  (0  <  A  <  1)  and 
■ J{Uk )  =  dF(Uk)/dU  the  Jacobian  matrix.  A  system  of  linear  block  tridiagonal  equations 
must  be  solved  at  each  iteration  for  corrections  to  the  previous  solution  vector.  In  the 
counterflow  diffusion  flame  problem,  the  cost  of  forming  (we  use  a  numerical  Jacobian)  and 
factoring  the  Jacobian  matrix  can  be  a  significant  part  of  the  cost  of  the  total  calculation. 
In  such  problems  we  apply  a  modified  Newton  method  in  which  the  Jacobian  is  re-evaluated 
only  periodically  [4]. 

The  solution  of  combustion  problems,  such  as  the  counterflow  diffusion  flame,  requires 
that  the  computational  mesh  be  determined  adaptively.  Many  of  the  methods  that  have 
been  used  to  determine  adaptive  grids  for  two-point  boundary  value  problems  can  be 
interpreted  in  terms  of  equidistributing  a  positive  weight  function  over  a  given  interval 
[5,6],  We  say  that  a  mesh  X  is  equidistributed  on  the  interval  [  —  L,  L)  with  respect  to  the 
non-negative  function  IV  and  the  constant  C  if 


dy  =  C, 


j  =  0, 1 . m  -  1. 


(2.57) 


We  determine  the  mesh  by  employing  a  weight  function  that  equidistributes  the  differ¬ 
ence  in  the  components  of  the  discrete  solution  and  its  gradient  between  adjacent  mesh 
points.  Upon  denoting  the  vector  of  N  =  K  - (-3  dependent  solution  components  by 
V  =  \U\.  U'i . Uy}T.  we  seek  a  mesh  X  such  that 


■y, 


dfj, 

dy 


dy 


< 


(U  max  U,  —  min  U, 

—  L<y<L  — L<y<L 


and 


j  =  0. 1 . rn  —  I 

i  =  1,2 . N 


d2U, 

~d\F 


dy 


< 


,  dU,  dU, 

•7  max  — - mm  — — 

dy  -i<v<l  dy 


j  —  1,2 . m  -  1 

i  =1.2 . .V 


(2.58) 


(2.59) 


where  d  and  */  are  small  numbers  less  than  one  and  the  maximum  and  minimum  values  of  U, 
and  dUJdy  are  obtained  from  a  converged  numerical  solution  on  a  previously  determined 

mesh. 


The  coarse  to  fine  grid  strategy  eliminates  many  of  the  convergence  difficulties  asso¬ 
ciated  with  solving  the  governing  equations  directly.  However,  convergence  of  Newton’s 
method  on  the  initial  grid  requires  a  “good”  initial  estimate  for  U°.  We  can  improve  the 
flame  sheet  starting  estimate  by  applying  a  time-dependent  starting  method  on  the  initial 
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grid.  We  remark  that,  fundamentally,  there  are  two  mathematical  approach*-  f. .r  .  1  v i n e 
one-dimensional  flame  problems  one  uses  a  transient  method  and  the  other  ->lv<-  tin 
steady-state  boundary  value  problem  directly.  Generally  speaking,  the  transient  meth 
ods  are  robust  but  computationally  inefficient  compared  to  the  boundary  value  method', 
which  are  efficient  but  have  less  desirable  convergence  properties.  Most  of  the  numerical 
techniques  that  have  been  used  to  solve  flame  problems  have  employed  a  time-dependent 
method.  Variations  of  this  approach  have  been  considered  by  a  variety  of  researchers  (see. 
e.g.,*{7-l l]).  In  these  methods,  the  original  nonlinear  two-point  boundary  value  problem  is 
converted  into  a  nonlinear  parabolic  mixed  initial-boundary  value  problem.  This  is  accom¬ 
plished  by  appending  the  term  ()(•)/ dt  to  the  left-hand  side  of  the  conservation  equations. 
This  results  in  a  semidiscrete  set  of  equations 


OU 

~dt 


=  F{U). 


(2. GO) 


with  appropriate  initial  conditions.  If  the  time  derivative  is  replaced,  for  example,  by  a 
backward  Euler  approximation,  the  governing  equations  can  be  written  in  the  form 

(Un+l  —  Un ) 

?(U',+l)  =  F(Un+l)  -  - — - l  =  0,  (2.61) 


where  for  a  function  g(t)  we  define  gn  =  g(tn )  and  where  the  time  step  r"+l  =  t"+'  - 
At  each  time  level  we  must  solve  a  system  of  nonlinear  equations  that  looks  very  similar 
to  the  nonlinear  equations  in  (2.55).  Newton’s  method  can  again  be  used  to  solve  this 
system.  The  important  difference  between  the  system  in  (2.55)  and  (2.61)  is  that  the 
diagonal  of  the  steady-state  Jacobian  is  weighted  by  the  quantity  l/rn+l.  This  produces  a 
better  conditioned  system  and  the  solution  from  the  time  step  ordinarily  provides  an 
excellent  starting  guess  to  the  solution  at  the  (n  +  l)8*'  time  level.  The  work  per  time  step 
is  similar  to  that  for  the  modified  Newton  iteration,  but  the  time-like  continuation  of  the 
numerical  solution  produces  an  iteration  strategy'  that  will,  in  general,  be  less  sensitive  to 
the  initial  starting  estimate  than  if  Newton’s  method  were  applied  to  (2.55)  directly.  As 
a  result,  when  we  ultimately  implement  Newton’s  method  on  the  steady-state  equations 
directly,  we  obtain  a  converged  numerical  solution  with  only  a  few  additional  iterations. 
This  time-dependent  starting  procedure  can  also  be  used  on  grids  other  than  the  initial 
one. 


2.4  Numerical  Results 

We  applied  the  flam*'  sheet  starting  estimate  to  a  diluted  methane-air  flame  in  the 
cylindrical  double-jet  configuration  (see  Figure  I).  As  we  discussed  in  Section  2.2,  the 
Hame  sheet  model  provides  initial  solution  profiles  for  the  mass  flux  in  the  transverse 
direction.  V.  the  similarity  function,  /'.  the  temperature.  T.  and  the  major  specie’s,  i.e. 
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CH.!,  (>2 ,  N2,  002  and  CO.  Gaussian  profiles  were  used  for  the  minor  species.  The  detailed 
kinetics  mechanism  used  in  the  calculations  is  listed  in  Table  1. 

After  getting  the  flame  sheet  starting  estimate,  we  solved  the  full  set  of  governing 
equations  in  a  tw'o-step  procedure.  Specifically,  we  determined  a  solution  to  the  mass, 
momentum,  and  species  equations  with  the  energy  equation  replaced  by  the  flame  sheet 
temperature  profile.  This  procedure  is  similar  to  the  two-pass  solution  method  used  in 
the  solution  of  adiabatic  premixed  laminar  flames  [3].  The  fixed  flame  sheet  temperature 
solution  (Tout)  was  then  used  to  obtain  a  solution  to  the  full  fluid  dynamic-thermochemistry 
model  (Tin).  This  procedure  helped  to  reduce  both  convergence  difficulties  and  total  CPU 
time. 

In  our  problem,  the  separation  distance  of  the  two  jets  was  1.4876  cm.  The  boundary 
conditions  at  the  fuel  jet  were  given  by 

V  =  2.8  x  10-2,  (2.62) 


(2.63) 

(2.64) 

(2.65) 

(2.66) 

(2.67) 

(2.68) 


The  mass  How  rate  was  in  units  of  gm/cm2-sec  and  the  densities  of  the  fuel  and  the  oxidizer 
mixtures  were  used  in  obtaining  the  value  of  the  similarity  function  at  the  fuel  jet.  The 
mass  flow  rate  boundary  condition  corresponds  to  a  fuel  duct  velocity  of  35  cm/sec.  The 
strain  rate  used  in  the  calculation  was  a  =  40  sec-1. 


A  solution  was  obtained  on  a  nonuniform  grid  consisting  of  38  grid  points.  This 
solution  was  then  used  as  the  starting  estimate  for  the  fixed  temperature  solution.  One 
hundred  adaptive  time  steps  were  taken  to  help  bring  the  solution  within  the  domain  of 
convergence  of  Newton  s  method  on  the  38  point  grid.  After  the  time  steps.  Newton’s 
method  converged  with  only  one  iteration.  Once  this  solution  was  obtained,  the  mesh  was 
refined  and  a  solution  was  calculated  on  the  finer  grid.  This  procedure  continued  until  the 
adaptive  mesh  criteria  were  satisfied.  The  refined  fixed  temperature  solution  was  then  used 
as  the  starting  estimate  for  the  complete  fluid  dynamic-thermochemistry  solution.  Two 
additional  grid  refinements  were  performed  to  obtain  a  final  solution  on  a  grid  consisting 
of  65  nonuniforin  points.  On  the  refined  grids  Newton's  method  converged  after  applying 
only  10-20  time  steps.  The  mesh  spacing  was  such  that  600  equispaced  points  would  have 
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been  needed  to  obtain  comparable  accuracy.  The  total  CPU  time  for  the  entire'  procedure 
was  approximately  three  hours  and  40  minutes  on  a  VAX-8000.  Approximately  0  CPU 
seconds  were  needed  for  the  flame  sheet  calculation,  130  minutes  for  the  Tout  calculation 
and  approximately  90  minutes  for  the  Tm  solution.  In  contrast,  we  were  unable  to  obtain  a 
complete  Tin  solution  for  this  problem  when  the  premixed  laminar  flame  starting  procedure 
was  employed. 

In  Figure  2  we  compare  the  flame  sheet  temperature  profile  with  the  calculated  finite 
rate  temperature  profile.  We  observe  that  the  agreement  is  generally  quite  good.  Similar 
results  are  illustrated  in  Figures  3  and  4  for  the  similarity  function  and  the  normal  velocity. 
The  flame  sheet  model  is  able  to  predict  adequately  the  location  of  the  stagnation  point 
along  with  the  "double  peak"  velocity  profile. 

With  the  success  of  the  flame  sheet  starting  procedure,  we  then  compared  the  results 
of  our  detailed  chemistry  model  with  the  experimental  measurements  of  Seshadri  [l].  In 
Figure  5  we  illustrate  the  experimental  (circles)  and  the  calculated  (solid  line)  temperature 
profiles  for  a  problem  with  similar  input  conditions  as  those  described  in  (2.62)-(2.68). 
We  observe  that  excellent  agreement  is  obtained  for  the  general  shape  of  the  profiles 
as  well  as  for  the  peak  temperatures  (maximum  experimental  =  1691  K  and  maximum 
computational  =  1680  I\).  On  the  rich  side  of  the  flame,  however,  the  temperature  profiles 
show  some  discrepancy.  These  differences  are  due  most  likely  to  the  neglect  of  C2  chemistry 
in  the  chemical  kinetics  mechanism.  In  Figures  6  and  7  we  compare  the  experimental 
and  computational  results  for  the  stable  major  chemical  compounds  with  the  physical 
coordinate  as  the  independent  variable.  The  solid  lines  represent  the  numerical  calculations 
and  the  points  represent  the  experimental  measurements.  The  profiles  reveal  excellent 
agreement  for  the  mole  fractions  of  CHA ,  02 ,  iV2 ,  H2 O  and  C02.  The  profiles  for  CO  and 
Hi.  however,  show  some  discrepancy.  The  differences  between  the  experiment  and  the 
calculations  are  due  again  to  the  neglect  of  C2  chemistry  in  the  calculations.  Minor  species 
and  radicals  are  illustrated  in  Figures  8  and  9. 

3.  FLAME  SHEET  MODEL  OF  AN  AXISYMMETRIC  DIFFUSION  FLAME 

The  modeling  of  axisymmetric  diffusion  flames  can  be  reduced  to  the  solution  of  a  set  of 
coupled  nonlinear  boundary  value  problems.  In  these  problems  we  desire  solution  profiles 
to  as  many  as  several  dozen  species  concentrations  in  addition  to  the  temperature  and  the 
velocity  fields.  Although  axisymmetric  flames  are  important  in  combustion  applications, 
they  have  received  relatively  little  attention  in  theoretical  flame  studies.  Part  of  this 
neglect  is  due  to  the  two-dimensional  nature  of  the  problem  coupled  with  the  complexities 
;i>M>ciated  with  the  combined  effects  of  transport  phenomena  and  chemical  processes.  In 
the  axisymmetric  diffusion  flame  we  consider,  a  fuel  jet  discharges  into  a  laminar  air  stream 
The  flames  ran  be  either  confined  or  unconfined.  In  both  cases  the  tubes  through  which 


the  fuel  and  the  oxidizer  flow  are  concentric  and  have  radii  R{  and  Ro,  respectively.  In 
both  configurations  the  two  gases  make  contact  at  the  outlet  of  the  inner  tube  and  a  flame 
resembling  a  candie  results. 


In  the  second  part  of  our  Phase  I  SBIR  we  investigated  the  generalization  of  our  one¬ 
dimensional  flame  sheet  model  to  two  dimensions.  The  goal  of  this  work  was  to  provide 
multi-dimensional  starting  estimates  for  the  full  chemistry-fluid  dynamic  model  that  will 
be  studied  in  Phase  II.  In  an  effort  that  complements  the  one  described  in  Section  4,  we 
used  a  stream  function-vorticity  as  opposed  to  a  primitive  variable  formulation.  In  this 
way  the  pressure  (and  the  continuity  equation  which  is  not  in  the  standard  elliptic  form 
of  the  remaining  equations)  could  be  eliminated  from  the  model  and  we  could  reduce  the 
size  of  the  system  to  be  solved. 


We  define  the  vorticity  as  the  amount  of  counter-clockwise  rotation  in  the  fluid.  We 


have 


dvr  dvz 


uvr  uuz 

“’-aT-aT'  (31) 

The  stream  function  t  is  used  to  replace  the  radial  and  axial  components  of  the  velocity 
vector  by  a  single  function.  It  is  defined  in  such  a  way  that  the  continuity  equation  is 
identically  satisfied.  We  have 

30 
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(3.2a) 

(3.26) 


Our  flame  sheet  model  follows  the  procedure  we  discussed  in  detail  in  Section  2.2.  We 
again  assume  a  single  global  irreversible  reaction.  The  diffusion  velocity  is  given  by  Fick’s 
law  and  the  heat  capacities  are  assumed  constant.  Our  model  can  be  written  in  the  form 


3/1  30  \  d  /  1  30  \ 
dz  \rp  dz  )  dr  \rp  dr  ) 


(3.3) 


r 


2 


'_3_ 

dz 


(  u>  d  il> 
\  r  dr 


d  /  u  dil>  \ 
dr  \r  dz  ) 


d 

dr 


i(H) 


dr 


v ?  +  v\ 


iso  p  —  0,  (3.4) 


d_ 

dr 


+ 


d_ 

dz 


rM 

dr 


d_ 

dr 


d 


dz 


—  rWFuFw  =  0, 


(3.5) 


12 


-i(*£K(*£KC 

-  “  rWoVou>  = 

d  (v  d^\  .  3  fYPdj>\  d  ( 
dr  \  P  dz  )  3j  \  3r  /  dr  \ 


d  /  n  dV'o 

*  (rpDo^r 


rW0u0w  =  0, 


dz  V  dr 


rpDp- 


(rpDp~^)  +  rWPuP™  =  °’ 

-S  (r-S)*E  ('-■?) -S  ("“-¥) 

-Jh?)- 


’  d 

(Td-±\ 

d 

frW\l 

dz 

V  dr) 

dr 

[Tdz)\ 

WpUphp  +  W0u0h0  —  WpUphp  . 

— r - te  =  0,  (3.9) 

cp 

where  ti>  is  again  the  rate  of  progress  of  the  reaction.  If  we  follow  the  argument  in  (2.30)- 
(2.31),  we  find  that  each  of  the  Shvab-Zeldovich  variables 


ZF  =  Yf  -  Yf^  +  ^T -T„), 

(3.10) 

Zo  =  Yo-Yo„  +  C’Z°l°(T-Tm), 

(3.11) 

Z^  =  Y'-Y^-qZZ{T~T-)- 

(3.12) 

ZN  =  Yn  —  yNoo, 

(3.13) 

satisfies  the  differential  equation 
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(3.14) 


One  can  show  that  all  of  the  Z*  are  proportional  to  each  other  and  to  a  conserved  scalar 
S  that  satisfies  an  equation  similar  in  form  to  (3.14).  Relations  for  the  temperature  and 
the  major  species  follow  exactly  as  in  Section  2.2. 

We  applied  the  two-dimensional  flame  sheet  model  to  an  unconfined  axisvmmetric 
methane-air  flame  with  an  inner  jet  (fuel)  radius  of  0.254  cm  and  an  outer  jet  (air)  radius 


V.  \  i  (  r, 


of  2.54  cm.  The  velocities  of  the  fuel  and  the  oxidizer  were  10  cm/sec  and  the  temperature 
of  the  two  gases  was  298  K.  Utilizing  a  generalization  of  the  solution  procedure  discussed 
in  Section  2.3,  we  obtained  solution  profiles  for  the  temperature,  major  species  and  the 
velocities  on  a  21  x  41  grid.  Results  of  the  calculations  are  illustrated  in  Figures  10-15. 

4.  PRIMITIVE  VARIABLE  METHODS  IN  TWO  DIMENSIONS 

£n  this  section  we  report  on  experimental  modifications  to  the  TEACH  [12]  code  for 
solving  the  conservation  equations  of  reacting  flows  in  primitive  variable  form.  We  ex¬ 
amined  primarily  the  CPU-intensive  SIMPLE  (Semi-Implicit  Method  for  Pressure-Linked 
Equations)  algorithm  [13]  for  handling  the  momentum-pressure  coupling,  which  lies  at 
the  heart  of  TEACH ,  and  to  which  TEACH  reduces  in  the  laminar,  incompressible,  single¬ 
component  case. 

TEACH  is  a  multidimensional  steady -state  compressible  flow  algorithm  based  on  a 
low-order  control-volume  discretization  of  the  primitive  variable  equations  in  conservation 
form  over  staggered,  orthogonal,  tensor- product  grids.  In  addition  to  the  continuity  and 
momentum  equations,  a  number  (in  principal  arbitrary)  of  conserved  scalar  equations  can 
be  accommodated  by  TEACH.  These  may  include  an  energy  equation,  multiple  specie  equa¬ 
tions,  and  the  empirical  k  -  c  turbulence  model,  which  supplements  the  laminar  equation 
system  with  two  extra  transport  equations  of  the  same  generic  form  while  modifying  the 
laminar  diffusion  terms  in  the  remainder  of  the  equation  set. 

A  large,  sparsely  coupled  set  of  nonlinear  algebraic  equations  results.  The  nonlineari¬ 
ties  arise  directly  from  the  advective  terms  and  the  source  terms,  and  indirectly  through 
temperature-,  pressure-,  and  composition-dependent  laminar  transport  properties  and 
thermodynamic  coefficients,  and  the  velocity-dependent  turbulent  transport  properties. 
A  multistage  variation  of  the  block  Gauss-Seidel  method  is  used  to  solve  this  nonlinear 
system.  The  outermost  stage  consists  of  cycling  between  a  Poisson-like  pressure  correction 
equation  (derived  from  a  truncated  substitution  of  the  discrete  momentum  equations  into 
the  discrete  continuity  equation,  which  is  thus  eliminated)  and  the  transport  equations 
for  all  of  the  other  unknown  fields.  Within  this  latter  block,  the  equations  are  relaxed 
cyclically,  field-by-field.  Within  the  sub-blocks  at  the  level  of  the  individual  fields,  the 
updates  are  computed  in  a  block-line  fashion  so  that  a  tridiagonal  matrix  is  the  largest 
implicit  aggregate  in  the  overall  calculation.  In  practice,  under-relaxation  of  the  updates 
is  necessary. 

Variations  of  the  TEACH  code  abound.  The  version  with  which  we  make  our  com¬ 
parisons  is  a  revision  obtained  from  one  of  the  original  authors,  as  described  in  [12],  Our 
principal  test  problem  for  this  section  is  the  two-dimensional,  incompressible,  axisymmet- 
ric,  nonreacting  flow  in  a  suddenly  expanded  laterally  heated  channel,  the  hydrodynamics 
of  which  are  described  in  [14],  This  is  the  test  problem  described  in  [12]  and  supplied  with 


the  official  machine-readable  copy  of  the  code,  being  bound  up  therewith  in  a  non-modular 
way. 


4.1  Advantages  and  Disadvantages  of  TEACH 

As  a  base  on  which  to  build  a  detailed-kinetics  reacting  flow  solver,  TEACH  has  several 
advantages,  which  motivated  us  to  improve  its  convergence  properties: 

•  The  pressure,  which  is  in  practice  one  of  the  fields  of  greatest  interest,  is  readily 
available  as  one  of  the  unknowns  of  the  problem,  whereas  the  pressure  is  eliminated  in 
a  stream  function-vorticity  formulation. 

•  The  TEACH  algorithm  generalizes  in  principle  to  three  dimensions,  the  space  in  which 
all  real  engineering  problems  lie. 

•  The  widely  used  k  —  t  turbulence  model  is  built  in. 

•  The  wide  user  base  of  TEACH-\ ike  codes  would  seem  to  portend  a  high  interest  in  the 
ultimate  production  version  of  the  reacting  flow  solver,  and  relative  ease  of  transfer- 
ability  to  the  interested  community. 

There  are  also  some  major  disadvantages: 

•  TEACH  has  many  parameters  (such  as  underrelaxation  factors  and  iteration  limits) 
that  are  difficult  to  tune  efficiently,  and  poor  diagnostics  by  which  to  determine  ‘con¬ 
vergence"  . 

•  Convergence  is  slow,  being  asymptotically  linear  at  best,  even  for  “easy"  problems 
(those  with  few  species,  the  conservation  equations  for  which  are  dominated  by  the 
linearly  implicit  part  of  the  source  term,  or  by  advection  or  diffusion). 

•  Convergence  is  unreliable  in  “difficult"  problems  (those  with  many  species,  the  con¬ 
servation  equations  for  which  are  dominated  by  parts  of  the  source  term  not  treated 
implicitly). 

•  There  are  no  convenient  data  structures  built  into  TEACH  to  assist  the  user  in  organiz¬ 
ing  the  large  volume  of  physical  data  (kinetic,  transport,  and  thermodynamic)  which 
must  be  supplied  in  realistic  detailed  modeling. 


Because  of  these  disadvantages,  much  of  our  work  in  preparation  for  detailed  modeling  j 

in  two  dimensions  has  been  carried  out  in  a  context  other  than  the  TEACH  code,  as  ) 

described  above.  However,  since  there  is  ultimate  interest  in  a  primitive  variables  staggered  I 

grid  code,  we  addressed  the  first  two  disadvantages  to  TEACH ,  as  described  in  the  balance 
of  this  section. 

! 

4.2  Modifications  of  the  TEACH  Code  -  Outer  Loop  » 

15  ' 


The  weakness  at  the  heart  of  the  TEACH  algorithm  is  in  the  solution  of  the  discrete 
nonlinear  system.  Linearization  by  decoupling  at  the  field-by-field  level  provides  only  a 
linear  convergence  rate  for  the  outer  iteration.  Furthermore,  the  degree  of  underrelaxation 
found  necessary  in  practice  causes  this  rate  to  be  slow.  The  most  natural  modification  was 
to  implement  a  Newton-like  algorithm  for  the  outer  loop. 

New  ton  methods  are  desirable  for  their  asymptotically  superlinear  convergence  rates. 
They  have  the  big  disadvantage  compared  to  simpler  iterative  schemes,  such  as  that  of 
plain  TEACH,  that  they  require  formation  and  factorization  of  a  Jacobian  matrix.  YVe 
examined  a  recently  developed  Jacobian-free  Newton-like  method  known  as  the  nonlinear 
generalized  minimum  residual  method  (NLGMR).  A  precursor  of  this  method  may  be 
found  in  [15].  More  recently,  very  impressive  results  have  been  obtained  with  NLGMR  in 
the  computational  fluid  dynamics  context  [  16] .  The  incorporation  of  NLGMR  into  TEACH 
can  be  accomplished  in  a  very  modular  fashion,  leaving  most  of  the  original  code  fully 
intact.  To  be  specific,  the  NLGMR  algorithm  is  inserted  between  the  control  driver  of 
TEACH  and  the  calls  to  the  field-by-field  solvers. 

Let  the  action  of  one  outer  iteration  of  the  TEACH  algorithm  be  denoted  by 

un+t  —  M(u„),  (4.1) 

where  the  vector  u  €  Rv  represents  the  discrete  unknowns  of  all  of  the  fields  (u.  r,  p, 
etc.)  and  M  is  the  nonlinear  mapping  that  produces  the  (n  +  l)st  iterate  from  the  nth.  In 
the  TEACH  code.  M  represents  one  pass  through  the  routines  CALCU.  CALCV,  CALCP.  etc.,  in 
their  proper  cyclic  order.  The  converged  solution  u*  is  the  root  of  the  system  of  nonlinear 
equations 

F(u )  =  u  —  M(u)  =  0.  (4.2) 

With  the  specification  of  the  initial  iterate,  «,>,  the  TEACH  algorithm  has  the  form  of  a 
nonlinear  block  Gauss-Seidel  method  for  F(u)  =  0. 

NLGMR  provides  a  means  of  accelerating  the  TEACH  iterations  as  follows  (see  [16]). 
Throughout  we  denote  by  ./(«)  the  Jacobian  matrix  of  F  evaluated  at  u.  When  there  is  no 
ambiguity,  ./(«)  will  be  denoted  simply  by  ./.  Suppose  that  un  is  the  current  approximation 
and  that  we  wish  to  find  a  new  approximation  of  the  form  u„+i  =  u„  +  b. 

We  write  b  in  the  form 

m 

f>  =  ft.-t’.a  (4.3) 

i  =  i 

where  the  v ,  £  Rv  are  ni  orthonormal  vectors,  the  specification  of  which  will  be  given 
shortly,  and  the  a,  are  unknowns  to  be  determined  at  each  step.  Ideally,  we  would  carry 
out  the  minimization  of 

l|F(Hn+*)|| 
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over  all  vectors  <■»  to  get  the  new  iterate  m„+i  =  u„  +  6.  Instead  of  solving  this  numerically 
difficult  problem  directly,  we  can  linearize  F  about  t/„  by  writing  F(un  +  d)  s  F(u„)  + 
./(«„)/>  and  seek  to  minimize 

||FK)  +  ^li 

over  all  vectors  d  of  the  form  (4.3).  If  we  choose  the  v,  to  form  an  orthonormal  basis  for  the 
Krylov  subspace  based  on  the  Jacobian,  namely  span  {t'i,  Jv\,  J2i\, . . . ,  Jm~l  t'i },  where 
»'i  =  F(Hr,)/||F(un)||.  then  this  minimization  can  be  carried  out  by  performing  tn  steps  of 
the  GMRES  algorithm  [17]  applied  to  the  linear  system 

Jb  =  —  F(u„),  (4.4) 

with  the  initial  iterate  b0  =  0.  Note  that  an  exact  solution  of  (4.4)  would  yield  the  Newton 
direction  .J~l  F{u„).  We  leave  most  of  the  details  of  GMRES  to  the  references;  however,  we 
note  several  of  its  important  properties.  If  rn  =  N  the  method,  though  iterative  in  nature, 
delivers  the  exact  solution  of  (4.4)  and  is  thus  a  strict  Newton’s  method  in  this  limit.  If 
the  Jacobian  is  suitably  conditioned,  however,  GMRES  usually  converges  for  m  <C  .V, 
resulting  in  considerable  savings  over  the  linear  algebra  expense  of  a  direct  elimination 
Newton  algorithm.  For  sufficiently  large  m,  its  storage  requirements  are  roughly  half 
those  of  theoretically  equivalent  Krylov  subspace  methods,  such  as  GCR  and  ORTHODIR 
[  18 j .  and  its  CPU  time  requirements  are  about  one-third  less.  Finally,  the  method  cannot 
break  down  before  finding  the  solution  (in  the  absence  of  roundoff  error)  for  arbitrary 
nonsy mmet ric  indefinite  matrices  J. 

Perhaps  the  most  important  aspect  of  the  NLGMR  algorithm  described  above  is  that 
the  Jacobian  matrix  .7  is  never  needed  explicitly.  The  only  operations  using  the  Jacobian 
matrix  .7  that  are  actually  required  in  the  implementation  are  matrix-vector  multiplications 
of  the  form  w  =  Jr,  which  can  be  approximated  by  finite  differences  of  F,  viz.. 

F(u  +  hv)  -  Fju) 

h  ( 

where  u  is  the  point  where  the  Jacobian  is  being  evaluated  and  h  is  some  carefully  chosen 
-mail  scalar.  In  comparison  to  the  typical  cost  of  evaluating  F(u),  which  requires  one 
application  of  the  mapping  A/(«).  there  is  very  little  overhead  in  the  NLGMR  algorithm, 
and  its  overall  cost  can  be  estimated  in  units  of  F-evaluations,  which  are  the  same  as  plain 
TEACH  outer  iterations. 

An  unfortunate  and  well-known  flaw  of  Newton's  method  is  that  the  domain  of  con¬ 
vergence  is  typically  quite  small.  In  practice,  NLGMR  invariably  diverges  when  applied 
to  TEACH  in  the  form  (4.2).  In  contrast,  the  TEACH  iterations  can  often  be  stabilized  by 
sufficient  underrelaxation  if  one  can  afford  the  price  of  the  concommitant  slow  convergence. 
A  reliable  cure  for  NLGMR.  following  [16],  is  to  introduce  damping  into  F(it).  that  is.  to 
replace  (4.2)  at  each  step  with 

F(u)  =  u  -  (1  -  A) A/(m)  -  AA/(u„),  (4.6) 
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where  A  is  a  (lamping  parameter  to  be  selected  adaptively,  0  <  A  <  1.  Selecting  A  =  0 
corresponds  to  Newton's  method.  Selecting  A  =  1  corresponds  to  ordinary,  unaccelerated 
TEACH.  If  A  is  held  at  1  until  the  TEACH  iteration  residuals  begin  to  exhibit  monotonic 
convergence,  a  subsequent  power-law  decrease  in  A  (as  a  function  of  iteration  count)  gen¬ 
erally  brings  about  asymptotically  superlinear  convergence.  To  be  specific,  let  p„  denote 
some  scaled  Euclidean  norm  of  the  components  of  ir(u„),  and  let  n  denote  the  iteration 
level  at  which  the  monotonic  residual  decrease  criterion  is  first  satisfied.  (In  practice, 
we  have  been  requiring  an  absolute  decrease  for  five  successive  iterations,  because  of  the 
notorious  oscillatory  convergence  profiles  of  plain  TEACH.)  Let  pn  denote  the  residual  at 
iteration  n.  Thereafter,  we  take 

A  n  =  (pn/pn)2-bn-\  (4.7) 

where  b  is  a  tuning  parameter  governing  the  rate  of  decrease  of  A,  supplied  by  the  user. 
Typically,  we  choose  b  in  the  range  0.7-0. 8,  which  promotes  a  fairly  rapid  transition  from 
plain  TEACH  to  Newton's  method  once  the  tail  of  the  TEACH  iterations  has  been  reached. 
The  pre-exponential  factor  is  designed  to  hasten  adaptively  the  transition  to  Newton’s 
method  according  to  the  rate  of  decrease  of  the  residual. 

In  Figures  16-19  we  show  the  type  of  savings  that  are  achievable  by  w'rapping  NLGMR 
around  TEACH  in  this  manner.  The  four  figures  are  in  two  pairs — one  pair  (Figures  16 
and  17)  for  a  16  x  16  grid  (14  interior  cells),  and  one  (Figures  18  and  19)  for  a  30  x  30  grid 
(28  interior  cells),  or  twice  the  resolution  of  the  first.  Each  pair  contrasts  the  convergence 
history  of  plain  and  NLGMR-accelerated  TEACH  for  a  laminar  version  of  the  standard  test 
problem  by  plotting  the  mass  and  momentum  equation  residuals  (solid  curves  -  individually 
nondimensionalized  as  in  the  original  TEACH)  and  the  composite  NLGMR  residual  (dashed 
curve)  against  the  number  of  evaluations  of  M{u).  In  the  16  x  16  case,  the  monotonically 
decreasing  tail  is  reached  at  33  function  evaluations.  Thereafter,  plain  TEACH  requires  105 
outer  iterations  (costing  one  function  evaluation  each)  to  reach  the  convergence  criterion 
of  a  relative  residual  reduction  of  10-8,  whereas  NLGMR-accelerated  TEACH  requires  only 
7  additional  outer  iterations  (costing  a  total  of  67  function  evaluations).  Each  plateau  in 
the  NLGMR  plots  represents  one  outer  iteration  during  which  several  GMRES  steps  may 
be  required.  In  the  30  x  30  case,  the  monotonically  decreasing  tail  is  reached  at  60  function 
evaluations.  Thereafter,  plain  TEACH  requires  153  outer  iterations  (costing  one  function 
evaluation  each)  to  reach  the  convergence  criterion  of  a  relative  residual  reduction  of  10-8. 
whereas  NLGMR-accelerated  TEA CH  requires  only  9  additional  outer  iterations  (costing  a 
total  of  108  function  evaluations).  In  terms  of  the  overall  number  of  function  evaluations 
from  start  to  finish  NLGMR  allows  reductions  of  28°o  and  21°o,  respectively.  In  no  case 
was  a  Krylov  subspace  of  dimension  higher  than  20  required. 

This  is  not  a  dramatic  improvement;  on  the  other  hand,  it  has  been  achieved  with 
virtually  no  extra  programming  cost  to  the  user,  and  some  possible  side  benefits  of  NLGMR 
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are  as  yet  unexplored  in  connection  with  TEACH.  For  instance,  in  [16]  the  authors  found 
it  possible  to  play  loosely  with  the  tuning  parameters  of  an  algorithm  being  accelerated 
by  NLGMR,  since  NLGMR  had  a  stabilizing  effect  on  the  overall  iteration.  Hence,  some 
of  the  uncertainty  that  accompanies  the  selection  of  TEACH  parameters  for  problems  of  a 
previously  unexplored  scale  may  be  relieved. 

We  conclude  that  the  damped  NLGMR  technique  provides  a  useful  enhancement  to 
plaiq.  TEACH  in  its  present  state  of  development.  Further  efforts  to  successfully  marry  these 
two  algorithms  (exogenous  to  our  AFOSR-sponsored  reacting  flow  work)  are  currently 
underway  by  SCA  consultants. 

4.3  Modifications  of  the  TEACH  Code  -  Inner  Loop 

In  view  of  the  preservation  of  the  TEACH  inner  loops  (as  embodied  in  the  routines 
CALCU.  CALCV.  CALCP,  etc.)  in  the  NLGMR-enhanced  version  described  above,  it  was  natural 
to  look  for  further  savings  in  solving  the  system  of  linear  equations  for  the  updates  to  each 
field.  Labeling  each  of  u,  v ,  p.  etc.  by  the  generic  field  0,  in  turn,  these  linear  systems  at 
each  inner  iteration  of  the  TEACH  algorithm  nominally  have  the  form 

ap4>p  =  ^o*0*  +  sp,  (4.8) 

* 

at  each  point  P  of  the  grid  on  which  0  is  defined,  where  the  sum  extends  over  the  four 
neighbors  of  P.  The  coefficients  arise  from  the  control-volume  discretization  of  the  advec- 
tion.  diffusion  and  source  operators  using  the  “best  available”  velocity  and  property  data. 
Because  the  advective  fluxes  are  of  hybrid  upwinded  form,  all  of  the  a*  are  nonnegative. 
Because  of  the  conservation  form  of  the  discretization, 

op  =  y (4.9) 
* 

In  fact,  the  field  updates  are  generally  underrelaxed  due  to  the  nonlinearities  hidden 
in  the  coefficients  of  (4.8).  Thus,  if  0*  represents  the  vector  of  unknowns  calculated  from 
(4  8).  where  the  coefficients  are  based  on  the  nth  iterate  d>",  then 

<2>n+l  =  +  (1  -  fi)on.  (4.10) 


where  //.  0  <  [i  <  1  i.s  an  underrelaxation  parameter,  which  may  be  set  differently  for 
each  field  (Typical  values  of  //  are  0.5  for  the  velocities  and  1.0  for  the  pressure,  with  tin1 
other  conserved  scalars  somewhere  in  between.)  If  the  degree  of  underrelaxation  desired 
in  (4.10)  is  known  a  priori,  it  may  be  built  into  the  coefficients  of  a  modified  version  of 
(4  8)  with  advantage.  In  fact,  the  TEACH  routine  block-line  solver  LISOLV  actually  solves 


0p0P  =  ^  0*0*  + 


(4.11) 


'  flflt  1 V  L.\  Lf  V  a". 


V  *  •  *  ’»  *.**  s’*  O*  -**  * 


10 


directly  for  0,,+1,  where  n'r  =  ap/ft  and  Sp  =  sp  +  (1  -  / i)a*P4>n .  This  is  algebraically 
equivalent  to  the  two-step  process  (4.8)-(4.9),  but  is  numerically  advantageous,  because 
it  increases  the  diagonal  dominance  of  the  linear  system,  thus  enhancing  the  convergence 
rate  of  the  block-line  iterative  solver. 

For  systems  of  appreciable  size,  GMRES  with  some  form  of  approximate  factorization 
preconditioner  is  among  the  best  methods  for  problems  of  the  form  (4.8).  Hence,  the  block¬ 
line  solver  LISOLV  was  replaced  with  several  differently  preconditioned  GMRES  methods 
from  PCGPAK  [19],  The  results  were  disappointing,  with  GMRES  typically  requiring  2.0 
to  2.5  times  as  much  CPU  time  as  the  block-line  solver  for  the  same  degree  of  residual 
reduction  on  a  16  x  16  grid  with  the  full  turbulent  test  problem  mechanism  in  place.  The 
conclusion  of  this  brief  exercise  is  that  for  the  amount  of  underrelaxation  already  built 
into  TEACH  to  stabilize  the  outer  nonlinear  iterations,  the  most  elementary  linear  solvers 
perform  most  efficiently  for  the  inner  iterations.  Therefore,  we  recommend  retaining  LISOLV 
in  conjunction  with  NLGMR-  TEACH. 

5.  CONCLUSIONS 

We  have  demonstrated  the  feasibility  of  combining  detailed  transport  phenomena  with 
complex  chemistry  in  the  solution  of  chemically  reacting  flows.  We  have  illustrated  the 
effectiveness  of  flame  sheet  initialization  procedures  and  a  time  integration-Newton  non¬ 
linear  equation  solver  in  the  modeling  of  a  counterflow  diffusion  flame.  We  have  also  made 
improvements  to  the  outer  iteration  of  the  primitive  variable  TEACH  code.  Our  goals  in  a 
Phase  II  SBIR  will  be  to  generalize  these  ideas  to  enable  the  efficient  solution  of  axisym- 
metric  (laminar  and  turbulent)  diffusion  flames.  Based  upon  our  Phase  I  results,  we  are 
confident  that  our  approach  will  be  successful. 
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A 

d 

E 

1. 

C774  +  M  = 

077.3  +  77  4  A/ 

1.00E417 

0.000 

86000. 

o 

077,  4  02  = 

C77:,  +  7702 

7.90E413 

0.000 

56000. 

3. 

077.,  4  H  = 

077.3  +  772 

2.20E404 

3.000 
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-4. 

OT74  4  O  = 

0773  +  077 

1.60E406 

2.360 
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5. 

077,  4  OH  = 
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1.60E+06 

2.100 
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G. 

CH20  4  OH 

-  7700  4  77,0 

7.53E412 

0.000 
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CH20  4  H  = 
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3.31E414 
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8. 

077*0  4-  A/ 
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3.31E41G 
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10. 
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11. 

7700  4  A/  = 

=  77  4  OO  4  A/ 

1.60E414 

0.000 
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12. 

HCO  4  H  = 
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4.00E413 

0.000 
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13. 

HCO  4  O  - 

077  4  00 

1.00E413 
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14. 
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15. 
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1G. 

CO  4  OH  = 

002  4  77 

1.51E407 

1.300 

-758. 

17. 

CO  4  o2  = 

CO 2  4  O 

1.60E413 

o.oon 

41000. 

18. 

CH3  4  O  ,  = 

07730  4  O 

7.00E412 

0.000 

25G52. 

19. 

C  H30  +  M 

=  C H20  4  77  4  M 

2.40E-r  13 

0.000 

28812. 

20. 

CH30  4  H  = 

=  C H20  4  772 

2.00E413 

0.000 

0. 

21. 

CHaO  4  OH  —  CH20  4  H20 

1.00E413 

0.000 

0. 
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CH:kO  4  O  - 

=  077,0  4  077 

1.00E413 

0.000 

0. 

23. 

C  HR)  4  Oa 

=  C H20  4  77 02 

C.30E410 

0.000 

2G00. 

21. 

CHt  4  02  = 

(' H20  4  OH 

5.2DE4 13 

0.000 

34574. 

25. 

G’//;,  4  0  — 

CH20  4  77 

6.80E4 13 

0.000 

0. 

2G. 

CH,  4  OH  - 

=  07720  4  772 

7.50E412 

0.000 

0. 

27. 

ho2  4  oo  = 

=  C02  4  0  77 

5.80E4 13 

0.000 

22934. 

28. 

H2  4  02  -  2 

OH 

1.70E4 13 

0.000 

47780. 

29. 

OH  4  H2  — 

7720  4  77 

1.17E409 

1.300 

362C. 

30. 

H  4  O  ,  =  O//  4  O 

5.13E41G 

-0.81G 

16507. 
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TABLE  I.  (continued) 

Reaction  Mechanism  Rate  Coefficients  In  The  Form  kf  =  AT0exp{- E0/ fiT). 
Units  are  moles,  cubic  centimeters,  seconds,  Kelvins  and  calories/mole. 


REACTION 

A 

0 

E 

31. 

O  +  H2  -  OH  +  H 

1.80E+10 

1.000 

8826. 

32. 

H  +  02  +  M  —  H02  +  Ma 

2.10E+18 

-1.000 

0. 

33. 

H  +  02  +  02  —  H02  +  02 

6.70E+19 

-1.420 

0. 

34. 

H  4-  02  +  iV 2  ^  H02  -f-  N2 

6.70E+19 

-1.420 

0. 

35. 

OH  +  H02  —  H20  +  02 

5.00E+13 

0.000 

1000. 

36. 

H  +  H02  =  2  OH 

2.50E+14 

0.000 

1900. 

37. 

O  +  H02  -  0-2  +  OH 

4.80E+13 

0.000 

1000. 

38. 

2  OH  —  O  +  H20 

6.00E+08 

1.300 

0. 

39. 

H2  + M  -  H  +  H  +  Mk 

2.23E+12 

0.500 

92600. 

40. 

02  +  A/  —  O  +  O  ■+■  M 

1.85E+11 

0.500 

95560. 

41. 

H  +  OH  +  M  —  H20  +  Mc 

7.50E+23 

-2.600 

0. 

42. 

H  +  H02  =  H2  +  02 

2.50E+13 

0.000 

700. 

43. 

ho2  +  ho2  —  h2o2  -t-  0 1 

2.00E+12 

0.000 

0. 

44. 

H202  +  M  -  OH  A  OH  +  M 

1.30E+17 

0.000 

45500. 

45. 

H202  +  H  =  H02  +  H2 

1.60E+12 

0.000 

3800. 

46. 

H202  +  OH  —  H20  +  H02 

1.00E+13 

0.000 

1800. 

Third  body  efficiencies:  k^(H20)  =  21k$(Ar),  k5(H2)  =  3.3Ar5(Ar),  k5{N2)  =  kr,(02)  =  0. 
Third  body  efficiencies:  kl2(H20)  =  6k\2{Ar),  kl2(H)  =  2kl2(Ar),  kl2(H2)  =  3ku(Ar). 
Third  body  efficiency:  kl4(H20)  =  20fci4(Ar). 
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Figure  2:  Comparison  between  the  calculated  flame  sheet 
temperature  profile  (dotted  line)  and  the  calculated  finite  rate 
chemistry  temperature  profile  (solid  line)  in  the  double-jet 
problem. 
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Figure  3:  Comparison  between  the  calculated  flame  sheet 
similarity  function  profile  (dotted  line)  and  the  calculated 
finite  rate  chemistry  similarity  function  profile  (solid  line)  in 
the  double-jet  problem. 


Figure  4:  Comparison  between  the  calculated  flame  sheet 
normal  velocity  profile  (dotted  line)  apd  the  calculated  fi¬ 
nite  rate  chemistry  normal  velocity  profile  (solid  line)  in  the 
double-jet  problem. 
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Figure  5:  Comparison  between  measured  (o)  and  calculated 
values  (solid  line)  of  the  temperature  profile. 
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Figure  7:  Comparison  between  measured  H^O  (□),  CO2  (°), 
Hi  (+)  and  CO  (o)  profiles  and  corresponding  calculated 
values  (solid  lines). 
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X  *  2.540  (IN  CM) 


Figure  11:  Flame  sheet  CH4  (mass  fraction)  distribution  in 
a  two-dimensional  axisymmetric  diffusion  flame. 
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X  *  2.540  (IN  CM) 


Figure  14:  Flame  sheet  CO2  (mass  fraction)  distribution 
a  two-dimensional  axisymmetric  diffusion  flame. 
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Figure  15:  Velocity  (vector)  distribution  in  a  two-dimensional 
axisymmetric  diffusion  flame. 
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Figure  16:  Residuals  plotted  against  number  of  function  eval¬ 
uations  (approximately  proportional  to  CPU  time)  for  Plain 
TEACH  for  the  test  problem  on  a  16  x  16  grid. 


Momentum,  Mass  and  Composite  Residuals 


LEGEND : 

M  -  Continuity  Equation  Residual 
U,V  -  Momentum  Equat i on  Residuals 
dashed  -  Composite  Residual 


Figure  17:  Residuals  plotted  against  number  of  function  eval¬ 
uations  for  NLGMR  TEACH  for  the  test  problem  on  a  16  x  16 
grid. 


41 


. I  i  i  mm] _ ULLii. 


Mcrrentim,  Mass  and  Composite  Residuals 


». 


Plain  TEACH  Convergence  (30  by  30) 
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Figure  18:  Residuals  plotted  against  number  of  function  eval¬ 
uations  for  TEACH  for  the  test  problem  on  a  30  x  30  grid. 


Mcmentixn.  Mass  and  Ccnposile  Residuals 


tn0  NLGMR  TEACH  Convergence  (30  by  30) 

1 1  1  1  |  1  1  1  '  |  i  i  i  i  |  i  i  i  i  |  i — i — 

;M'l 

10"  I"  jJ'v. 

:  tf'V",. 


»0  100  150  2C 

Function  Evaluations 
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Figure  19:  Residuals  plotted  against  number  of  function  eval¬ 
uations  for  NLGMR  TEACH  for  the  test  problem  on  a  30x  30 
grid. 
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